function [wave,period,scale,coi,power,sig95,global_ws,global_signif]=wavetest(dt,t0,pad,dj,j1,lag1,mother,sst)

%WAVETEST Example Matlab script for WAVELET, using NINO3 SST dataset
%
% See "http://paos.colorado.edu/research/wavelets/"
% Written January 1998 by C. Torrence
%
% Modified Oct 1999, changed Global Wavelet Spectrum (GWS) to be sideways,
%   changed all "log" to "log2", changed logarithmic axis on GWS to
%   a normal axis.
%
%
% ------------------------------------------------------------------------
% Modified Oct 2006 (G. Charria)
%
%
%
% INPUT: dt: step (in time or space) 
%
%        t0: value of the first time step or the first position
%
%        pad: 1 - pad the time series with zeroes (recommended)
% 
%        dj: this will do 1/dj sub-octaves per octave
%
%        j1: xx/dj, this says do xx powers-of-two with dj sub-octaves each
%
%        lag1: lag-1 autocorrelation for red noise background
%
%        mother: type of wavelet 
%        - a string, equal to 'MORLET' or 'PAUL' or 'DOG'
%
%        sst: data to analyse
%
% OUTPUT: 
%        wave: complex wavelet coefficient
%
%        period: period for a temporal analysis or wavelength for a
%                spatial analysis. The vector of "Fourier" periods 
%                (in time units) that corresponds to the Scales.
%        
%        scale: the vector of scale indices, given by S0*2^(j*DJ), 
%               j=0...J1 where J1+1 is the total # of scales.
%
%        coi: The Cone-of-Influence, which is a vector of N points 
%             that contains the maximum period of useful information
%             at that particular time.
%             Periods greater than this are subject to edge effects.
%
%        power: (abs(wave)).^2 - wavelet power spectrum
%
%        sig95: where ratio > 1, power is significant
%
% -----------------------------------------------------------------




%clear all
%close all

%load 'sst_nino3.dat'   % input SST time series
%sst = sst_nino3;

%load 'sla.dat'
%sst=sla(:,2);


%sst=1e15*sin((1:504)*pi);
%sst=sin((1:504))+4.;
ymean=mean(sst);
old_sst=sst;


%------------------------------------------------------ Computation
%dt = 0.25 ;
%pad = 1;      % pad the time series with zeroes (recommended)
%dj = 0.25;    % this will do 4 sub-octaves per octave
%j1 = 9/dj;    % this says do 7 powers-of-two with dj sub-octaves each
%lag1 = 0.72;  % lag-1 autocorrelation for red noise background
%mother = 'Morlet';
% ----
%jdeb=1;
%jfin=j1;
% ----

variance = std(sst)^2;

n = length(sst);
% time = (0:length(sst)-1)*dt + t0 ;  % construct time array
% xlim = [time(1)-dt,time(n)+dt];  % plotting range
s0 = 2*dt;    % this says start at a scale of 6 months

% Wavelet transform:
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (abs(wave)).^2 ;        % compute wavelet power spectrum

% Significance levels: (variance=1 for the normalized SST)
[signif,fft_theor] = wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother);
sig95 = (signif')*(ones(1,n));  % expand signif --> (J+1)x(N) array
sig95 = power ./ sig95;         % where ratio > 1, power is significant

% Global wavelet spectrum & significance levels:
global_ws = variance*(sum(power')/n);   % time-average over all times
dof = n - scale;  % the -scale corrects for padding at edges
global_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);

% Scale-average between El Nino periods of 2--8 years
%avg = find((scale >= 2) & (scale < 8));
Cdelta = 0.776;   % this is for the MORLET wavelet
%scale_avg = (scale')*(ones(1,n));  % expand scale --> (J+1)x(N) array
%scale_avg = power ./ scale_avg;   % [Eqn(24)]
%scale_avg = variance*dj*dt/Cdelta*sum(scale_avg(avg,:));   % [Eqn(24)]
%scaleavg_signif = wave_signif(variance,dt,scale,2,lag1,-1,[2,7.9],mother);

%whos

%------------------------------------------------------ Plotting
% figure;
% %--- Plot time series
% subplot('position',[0.1 0.75 0.65 0.2])
% plot(time,sst)
% set(gca,'XLim',xlim(:))
% xlabel('Time (year)')
% ylabel('NINO3 SST (degC)')
% title('a) NINO3 Sea Surface Temperature (seasonal)')
% hold off

% %--- Contour plot wavelet power spectrum
% subplot('position',[0.1 0.37 0.65 0.28])
% levels = [0.0625,0.125,0.25,0.5,1,2,4,8,16] ;
% Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
% %contour(time,log2(period),log2(power),log2(levels));  %*** or use 'contourfill'
% imagesc(time,log2(period),log2(power));  %*** uncomment for 'image' plot
% xlabel('Time (year)')
% ylabel('Period (years)')
% title('b) NINO3 SST Wavelet Power Spectrum')
% set(gca,'XLim',xlim(:))
% set(gca,'YLim',log2([min(period),max(period)]), ...
% 	'YDir','reverse', ...
% 	'YTick',log2(Yticks(:)), ...
% 	'YTickLabel',Yticks)
% % 95% significance contour, levels at -99 (fake) and 1 (95% signif)
% hold on
% contour(time,log2(period),sig95,[-99,1],'k');
% hold on
% % cone-of-influence, anything "below" is dubious
% plot(time,log2(coi),'k')
% hold off

% %--- Plot global wavelet spectrum
% subplot('position',[0.77 0.37 0.2 0.28])
% plot(global_ws,log2(period))
% hold on
% plot(global_signif,log2(period),'--')
% hold off
% xlabel('Power (degC^2)')
% title('c) Global Wavelet Spectrum')
% set(gca,'YLim',log2([min(period),max(period)]), ...
% 	'YDir','reverse', ...
% 	'YTick',log2(Yticks(:)), ...
% 	'YTickLabel','')
% set(gca,'XLim',[0,1.25*max(global_ws)])

% % %--- Plot 2--8 yr scale-average time series
% % subplot('position',[0.1 0.07 0.65 0.2])
% % plot(time,scale_avg)
% % set(gca,'XLim',xlim(:))
% % xlabel('Time (year)')
% % ylabel('Avg variance (degC^2)')
% % title('d) 2-8 yr Scale-average Time Series')
% % hold on
% % plot(xlim,scaleavg_signif+[0,0],'--')
% % hold off

% % end of code




